/* 
 * File:   tools.h
 * Author: reinhard
 *
 * Created on 30. März 2011, 13:30
 */

#ifndef TOOLS_H
#define	TOOLS_H

#ifdef	__cplusplus
extern "C" {
#endif
#include "matrix.h"
#include "poly.h"
#include "remaindering.h"
#define MAXN 15
#define MYMOD(a,b) ((a)%(b)+(b))%(b)
#include <gmp.h>
#include <stdlib.h>
    int n;
    int adj[MAXN][MAXN];
    int pow2; //2 to the power of n
    int *p; //p value from bhkk-paper
    int *c; //c value from bhkk-paper
    int** P_U; //polynomial according to paper of björklund et.al. : size at first position, exponents in increasing order
    int** C_U;
    int* P_U_I; //[MAXN*MAXN + 1]; //P-polynomial to the power of i
    int* C_U_J; //[MAXN*MAXN + 1]; //C-polynomial to the power of J
    int* PC_IJ; //[2 * MAXN*MAXN*MAXN + 1]; //(P^i)*(C*j)
    int C_D[MAXN + 1][MAXN + 1]; //actual coefficients
    int* w_flat;
    int* new_w_flat;
    int* a_flat;
    int* b_flat;
    int* c_flat;
    int module;
    // definitions for more dimensional arrays

#define NewW(i,j,l) new_w_flat[(i)*n*n+(j)*n+l]
#define A(S,j) a_flat[(S)*(n+1) + (j)]
#define B(S,j) b_flat[(S)*(n+1) + (j)]
#define ELEM(i,S) ((1 << (i)) & S)
#define C_FLAT(S,j,i) c_flat[(S)*(n+1)+(j)][(i)]
#define C_POLY(S,j)   c_flat[(S)*(n+1)+(j)]	
#define myPow(x) (1 << (x))

    typedef int* polynom;
    typedef int set;

    void parseGraph(char* filename) { //routine to parse graph from file (specified format in wiki)
        FILE* f = fopen(filename, "r");
        if (f == NULL) {
            fprintf(stderr, "could not open file %s\n", filename);
            exit(EXIT_FAILURE);
        }
        int e;
        int i = 0;
        int j = 0;
        fscanf(f, "%d", &n);
        if (n > MAXN)fprintf(stderr, "exceeded MAXN, please adjust MAXN in tools.h\n");
        while (!feof(f)) {
            fscanf(f, "%d", &e);
            if (!isspace(e)) {
                if (e < 0 && e > 1) {
                    fprintf(stderr, "value different from 1 or 0 in adjacency matrix\n");
                    exit(EXIT_FAILURE);
                } else {
                    adj[j][i] = e;
                    j = j + (i == (n - 1) ? 1 : 0);
                    i = (i + 1) % n;
                }
            }
        }
        fclose(f);
    }

    void freeAll() { //cascade for cleaning up Memory   
        int S, i, j, l;
        for (S = 0; S < pow2; S++) {
            free(C_U[S]);
            free(P_U[S]);
        }
        free(C_U);
        free(P_U);
        free(c);
        free(p);
        free(PC_IJ);
        free(P_U_I);
        free(C_U_J);
    }

    void init() {
        P_U_I = (int*) malloc(((n + 1)*(n + 1) + 1) * sizeof (int));
        P_U_I[0] = 0;
        C_U_J = (int*) malloc(((n + 1)*(n + 1) + 1) * sizeof (int));
        C_U_J[0] = 0;
        PC_IJ = (int*) malloc(((n + 1)*(n + 1)*(n + 1) + 1) * sizeof (int));
        PC_IJ[0] = 0;
    }

    int getPos(int node, int S) { //get Position of node in subset S if reduced to induced graph
        if (S == 0)return 0;
        int i = 0;
        int count = 0;
        int mask = 1;
        while (mask != (1 << node)) {
            if (mask & S)
                count++;
            mask = mask << 1;
        }
        //	printf("counted position for node %d in set of size %d: %d \n",node,count, getNumOfNodes(S));
        return count;
    }

    void polyCopy(polynom a, polynom r) { //copy polynomial a to memory location r
        int i;
        for (i = 0; i <= a[0]; i++)
            r[i] = a[i];
    }

    void getSubsets(int *result, int superSet, int n) { // avoid iterating over non-subsets of the superset by shrinking, iterating over all subsets and memorizing positions of elements
        //printf("computing subsets for %d\n",set);
        int i = 0, j = 0;
        int shrunkSet = myPow(n);
        int* elements = (int*) malloc(n * sizeof (int));
        while (superSet > 0) {
            if (1 & superSet) {
                elements[j++] = i;
            }
            superSet = superSet >> 1;
            i++;
        }
        set S;
        for (S = 0; S < shrunkSet; S++) {
            result[S] = expandSubset(S, n, elements);
            //	printf("computed Subset %d\n",result[S]);
        }
        free(elements);
        //	printf("done computing subsets for %d\n",set);
    }

    int expandSubset(int subSet, int n, int* positions) { // from shrunk representation and positions, reconstruct corresponding subsets produced by getSubsets()   
        int set = 0;
        int i;
        for (i = 0; i < n; i++) {
            set = set | (((1 << i) & subSet ? 1 : 0) << positions[i]);
        }
        return set;
    }

    void getSubMatrix(int nodes, int* result) { //extract submatrix from given matrix to result
        int i, j, k, l;
        k = 0;
        l = 0;
        int size = getNumOfNodes(nodes);
        for (i = 0; i < n; i++) {
            if (ELEM(i, nodes)) {
                for (j = 0; j < n; j++) {
                    if (ELEM(j, nodes)) {
                        matrix(result, k, l, size) = adj[i][j];
                        l++;
                        //				printf("l : %d, size: %d \n",l, size);												
                    }
                }
                k++;
                //	printf("k : %d \n",k);
            }
            l = 0;
        }
    }

    int getNumOfNodes(set U) { //get size of vertex subset by adding number of 1s in binary representation
        int n = 0;
        int S = U;
        while (S > 0) {
            n += (S & 1);
            S = S >> 1;
        }
        return n;
    }

    int popCnt(set U) {
        switch (sizeof (set)) {
            case 4:
                U = U - ((U >> 1) & 0x55555555);
                U = (U & 0x33333333) + ((U >> 2) & 0x33333333);
                U = (U & 0x0F0F0F0F) + ((U >> 4) & 0x0F0F0F0F);
                U = (U & 0x00FF00FF) + ((U >> 8) & 0x00FF00FF);
                return (U & 0x0000FFFF) + ((U >> 16) & 0x0000FFFF);
            default:
                return getNumOfNodes(U);
        }
    }

    int mypow(int x) { //return 2 to the power of x
        return (1 << x);
    }

    int fac(int x) { //straight forward computation of faculty
        int i;
        int fac = 1;
        for (i = 1; i <= x; i++)
            fac *= i;
        return fac;
    }

    int binCoeff(int n,int k){
	//printf("%d choose %d = %d\n",n,k, fac(n)/(fac(k)*fac(n-k)));
	return fac(n)/(fac(k)*fac(n-k));
    }


#ifdef	__cplusplus
}
#endif

#endif	/* TOOLS_H */

